#include <life/options.hpp>
#include <life/lifealg/backend.hpp>
#include <life/lifediscr/functionspace.hpp>
#include <life/lifediscr/region.hpp>
#include <life/lifepoly/im.hpp>
#include <life/lifefilters/gmsh.hpp>
#include <life/lifefilters/exporter.hpp>
#include <life/lifepoly/polynomialset.hpp>
#include <life/lifevf/vf.hpp>

#include <life/lifefilters/geotool.hpp>

/** use Life namespace */
using namespace Life;
using namespace Life::vf;


inline
po::options_description
makeOptions()
{
  po::options_description code1doptions("Code1D options");
  code1doptions.add_options()
    ("hsize", po::value<double>()->default_value( 0.1 ), "mesh size")
    ("hsize1", po::value<double>()->default_value( 0.1 ), "mesh size1")
    ("hsize2", po::value<double>()->default_value( 0.1 ), "mesh size2")
    ("nu", po::value<double>()->default_value( 0.5 ), "nu")
    ("alphaX", po::value<double>()->default_value( 0.01 ), "alphaX")
    ("alphaY", po::value<double>()->default_value( -2.0 ), "alphaY")
    ("bccoeff", Life::po::value<double>()->default_value( 100.0 ), "coeff for weak Dirichlet conditions")
    ("toler", Life::po::value<double>()->default_value( 1e-6 ), "error tol")
    ("xc", Life::po::value<double>()->default_value( 0.0 ), "param de recouvrement des maillages")
    ("penaldir", Life::po::value<double>()->default_value( 10 ),
     "penalisation parameter for the weak boundary Dirichlet formulation")
    ("gmsh", po::value<float>()->default_value( 2.1 ), " version of gmsh( 2.0 or 2.1)")
    ;
  return code1doptions.add( Life::life_options() );
}

/**
 * This routine defines some information about the application like
 * authors, version, or name of the application. The data returned is
 * typically used as an argument of a Life::Application subclass.
 *
 * \return some data about the application.
 */
inline
AboutData
makeAbout()
{
  AboutData about( "MortarMethod" ,
		   "MortarMethod" ,
		   "0.2",
		   "Code MortarMethod",
		   Life::AboutData::License_GPL,
		   "Copyright (c) 2010 Universite Joseph Fourier");

  about.addAuthor("Vincent Chabannes", "developer", "vincent.chabannes@imag.fr", "");
  return about;

}


template<uint Dim,uint N,uint OrderGeo>
class Localisation
  :
  public Application
{
  typedef Application super;
public:

  //! Polynomial order \f$P_2\f$
  //    static const uint16_type Order = 5;
  ///static const uint Order = N;

  typedef double value_type;

  //---------------------------------------------------------------------------------------//

  typedef Backend<value_type> backend_type;
  typedef boost::shared_ptr<backend_type> backend_ptrtype;

  typedef backend_type::sparse_matrix_type sparse_matrix_type;
  typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
  typedef backend_type::vector_type vector_type;
  typedef backend_type::vector_ptrtype vector_ptrtype;

  //---------------------------------------------------------------------------------------//

  typedef Simplex<Dim,OrderGeo,Dim> convex_type;
  typedef Mesh<convex_type> mesh_type;
  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;

  //---------------------------------------------------------------------------------------//

  typedef Lagrange<N,Scalar,PointSetFekete> basis_u_type;
  typedef Lagrange<0,Scalar> basis_l_type;
  typedef bases<basis_u_type, basis_l_type> basis_type;

  typedef FunctionSpace<mesh_type, basis_type> space_type;
  typedef boost::shared_ptr<space_type> space_ptrtype;

  typedef typename space_type::element_type element_type;
  typedef typename element_type::template sub_element<0>::type element_0_type;
  typedef typename element_type::template sub_element<1>::type element_1_type;

  //---------------------------------------------------------------------------------------//

  typedef Exporter<mesh_type> export_type;
  typedef boost::shared_ptr<export_type> export_ptrtype;

  /**
   * Constructor
   */
  Localisation( int argc, char** argv, AboutData const& ad,
		po::options_description const& od )
    :
    super( argc, argv, ad, od ),
    M_backend( backend_type::build( this->vm() ) ),
    M_backend2( backend_type::build( this->vm() ) ),
    //M_backend( backend_type::build( BACKEND_TRILINOS ) ),
    meshSize( this->vm()["hsize"].template as<double>() ),
    meshSize1( this->vm()["hsize1"].template as<double>() ),
    meshSize2( this->vm()["hsize2"].template as<double>() ),
    exporter( Exporter<mesh_type>::New( this->vm(), this->about().appName() ) ),
    exporter1( Exporter<mesh_type>::New( this->vm(), "domain_1" ) ),
    exporter2( Exporter<mesh_type>::New( this->vm(), "domain_2" ) )
  {
  }

  mesh_ptrtype createMesh( std::string __name,std::string __str) ;

  mesh_ptrtype createMesh( double meshSize );
  mesh_ptrtype createMesh1( double meshSize ,double x_couvr);
  mesh_ptrtype createMesh2( double meshSize, double x_couvr);

    
  // run the application
  void run();
  void runMultiDomain();
  void runTest();
private:

  // solve the system \f$D u = F\f$
  void solve( sparse_matrix_ptrtype const& D, element_type& U, vector_ptrtype const& F, bool is_sym );
  //void solve( sparse_matrix_ptrtype& D, element_type& u, vector_ptrtype& F );
  void problem(element_0_type& u , element_0_type& v , sparse_matrix_ptrtype& A , vector_ptrtype& B );
  void problem1(element_type& U1 , element_type& U2  );
  void problem2(element_type& U2 , element_type& U1  );

  // export results to ensight format (enabled by  --export cmd line options)
  void exportResults( element_0_type& u );
  void exportResults( element_type& U1,element_type& U2,double time );

private:

  // linear algebra backend
  backend_ptrtype M_backend,M_backend2;
  sparse_matrix_ptrtype A1;
  vector_ptrtype B1;
  sparse_matrix_ptrtype A2;
  vector_ptrtype B2;
  // mesh characteristic size
  double meshSize,meshSize1,meshSize2;

  // exporter factory
  export_ptrtype exporter;
  export_ptrtype exporter1;
  export_ptrtype exporter2;
}; // Interpolate



template<uint Dim,uint N,uint OrderGeo>
typename Localisation<Dim,N,OrderGeo>::mesh_ptrtype
Localisation<Dim,N,OrderGeo>::createMesh( std::string __name,std::string __str)
{
  //typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
  //typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
  mesh_ptrtype mesh( new mesh_type );

  Gmsh gmsh;
  gmsh.setOrder(OrderGeo);

  std::string fname = gmsh.generate( __name, __str  );

  ImporterGmsh<mesh_type> import( fname );

  float gmsh_version = 2.1;//this->vm()["gmsh"].template as<float>();
  if (gmsh_version==2.0) import.setVersion( "2.0" );
  else if (gmsh_version==2.1) import.setVersion( "2.1" );

  mesh->accept( import );
  mesh->components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
  mesh->updateForUse();

  return mesh;

}


template<uint Dim,uint Order,uint OrderGeo>
typename Localisation<Dim,Order,OrderGeo>::mesh_ptrtype
Localisation<Dim,Order,OrderGeo>::createMesh1( double meshSize,double x_couvr )
{
  mesh_ptrtype mesh( new mesh_type );
  GeoTool::Node x1(-1,-1);
  GeoTool::Node x2(0,1);
  GeoTool::Rectangle R(meshSize,"Omega1",x1,x2);
  R.setMarker(_type="line",_name="Bord",_marker1=true,_marker3=true,_marker4=true);
  R.setMarker(_type="line",_name="Interface",_marker2=true);
  R.setMarker(_type="surface",_name="Omega1",_marker1=true);
      
  return createMesh("Omega1",R.geoStr());

}
//----------------------------------------------------------------------------------------------------------//
//----------------------------------------------------------------------------------------------------------//
//----------------------------------------------------------------------------------------------------------//

template<uint Dim,uint Order,uint OrderGeo>
typename Localisation<Dim,Order,OrderGeo>::mesh_ptrtype
Localisation<Dim,Order,OrderGeo>::createMesh2( double meshSize, double x_couvr )
{
  mesh_ptrtype mesh( new mesh_type );
  GeoTool::Node x1(0,-1);
  GeoTool::Node x2(1,1);
  GeoTool::Rectangle R(meshSize,"Omega2",x1,x2);
  R.setMarker(_type="line",_name="Bord",_marker1=true,_marker3=true,_marker2=true);
  R.setMarker(_type="line",_name="Interface",_marker4=true);
  R.setMarker(_type="surface",_name="Omega2",_marker1=true);
      
  return createMesh("Omega2",R.geoStr());

}


template<uint Dim,uint N,uint OrderGeo>
void
Localisation<Dim,N,OrderGeo>::exportResults( element_0_type& U )
{
  Log() << "exportResults starts\n";
  exporter->step(0)->setMesh( U.functionSpace()->mesh() );
  exporter->step(0)->add( "u", U );
  exporter->save();

} // Interpolate::export

template<uint Dim,uint Order,uint OrderGeo>
void
Localisation<Dim,Order,OrderGeo>::exportResults( element_type& U1, element_type& U2, double time  )
{
  
  element_0_type u1 = U1.template element<0>();
  element_1_type p1 = U1.template element<1>();

  element_0_type u2 = U2.template element<0>();
  element_1_type p2 = U2.template element<1>();

  exporter1->step( time )->setMesh( u1.functionSpace()->mesh() );
  exporter1->step( time )->add( "u1", u1 );
  exporter1->step( time )->add( "lambda1", p1 );
  exporter1->save();

  exporter2->step( time )->setMesh( u2.functionSpace()->mesh() );
  exporter2->step( time )->add( "u2", u2 );
  exporter2->step( time )->add( "lambda2", p2 );
  exporter2->save();

} // Test_Schwarz::exportResults

template<uint Dim,uint N,uint OrderGeo>
void
Localisation<Dim,N,OrderGeo>::run()
{

  if ( this->vm().count( "help" ) )
    {
      std::cout << this->optionsDescription() << "\n";
      return;
    }

  using namespace Life::vf;  
 
  //-----------------------------------------------------------------------------------//
#if 0
  int MyOrder=N;

  this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
			  % this->about().appName()
			  % convex_type::name()
			  % MyOrder
			  % this->vm()["hsize"].template as<double>()
			  );

  //-----------------------------------------------------------------------------------//

  mesh_ptrtype mesh = createMesh( meshSize );

  space_ptrtype Xh = space_type::New( mesh );

  element_type u( Xh, "u" );
  element_type v( Xh, "v" );

  //-----------------------------------------------------------------------------------//

  AUTO( SolAnalytique, cos(M_PI*Px()/2.0)*cos(M_PI*Py()/2.0) );

  vector_ptrtype R( M_backend->newVector( Xh ) );
  sparse_matrix_ptrtype J( M_backend->newMatrix(Xh,Xh) );

  this->problem(u,v,J,R);

  exportResults(u);
#endif
}

template<uint Dim,uint Order,uint OrderGeo>
void
Localisation<Dim,Order,OrderGeo>::runMultiDomain()
{

  if ( this->vm().count( "help" ) )
    {
      std::cout << this->optionsDescription() << "\n";
      return;
    }

  using namespace Life::vf;  
 
  //-----------------------------------------------------------------------------------//

  //int MyOrder=N;

  this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
			  % this->about().appName()
			  % convex_type::name()
			  % Order
			  % this->vm()["hsize"].template as<double>()
			  );

  //-----------------------------------------------------------------------------------//

  double xc = this->vm()["xc"].template as<double>();
  mesh_ptrtype mesh1 = createMesh1( meshSize1,xc );
  mesh_ptrtype mesh2 = createMesh2( meshSize2,xc );


  space_ptrtype Xh1 = space_type::New( mesh1 );
  space_ptrtype Xh2 = space_type::New( mesh2 );

  element_type U1( Xh1, "U1" );
  element_0_type u1 = U1.template element<0>();
  element_1_type p1 = U1.template element<1>();

  element_type U2( Xh2, "U2" );
  element_0_type u2 = U2.template element<0>();
  element_1_type p2 = U2.template element<1>();

  element_type U1old( Xh1, "U1old" );
  element_0_type u1old = U1old.template element<0>();
  element_1_type p1old = U1old.template element<1>();

  element_type U2old( Xh2, "U2old" );
  element_0_type u2old = U2old.template element<0>();
  element_1_type p2old = U2old.template element<1>();

  //double __ppp =math::sqrt( integrate( elements(mesh1), _Q<2*Order>(), idv(u1)*idv(u1) ).evaluate()(0,0) );
  //double __ppp2 =math::sqrt( integrate( elements(mesh2), _Q<2*Order>(), idv(u2)*idv(u2) ).evaluate()(0,0) );
  
  //-----------------------------------------------------------------------------------//

  A1=M_backend->newMatrix(Xh1,Xh1);
  B1=M_backend->newVector( Xh1 );

  A2= M_backend->newMatrix(Xh2,Xh2);
  B2= M_backend->newVector( Xh2 );

  //-----------------------------------------------------------------------------------//

  AUTO( SolAnalytique, cos(M_PI*Px()/2.0)*cos(M_PI*Py()/2.0) );
  
  double erreur=1;
  double err_tol=this->vm()["toler"].template as<double>();
  double err1,err2,norm1,norm2;
  int cpt=0;
  std::cout <<"\nola\n";    
  while (/*erreur>=err_tol&&*/ cpt <10) 
    {
   
      //U1old=U1;
      u1old=u1;
      p1old=p1;

      this->problem1(U1,U2);
      norm1 =math::sqrt( integrate( elements(mesh1), _Q<2*Order>(), idv(u1)*idv(u1) ).evaluate()(0,0) );
      err1=math::sqrt( integrate(elements(mesh1), _Q<2*Order>(), (idv(u1)-idv(u1old))*(idv(u1)-idv(u1old)) ).evaluate()(0,0) );
      if (norm1!=0)  err1=err1/norm1;

      //U2old=U2;
      u2old=u2;
      p2old=p2;
      this->problem2(U2,U1);
      norm2 =math::sqrt( integrate( elements(mesh2), _Q<2*Order>(), idv(u2)*idv(u2) ).evaluate()(0,0) );
      err2=math::sqrt( integrate(elements(mesh2), _Q<2*Order>(), (idv(u2)-idv(u2old))*(idv(u2)-idv(u2old)) ).evaluate()(0,0) );
      if (norm2!=0)  err2=err2/norm2;

      erreur=err1+err2;
      std::cout<<"\n erreur = "<<erreur<<" ; err1 = "<<err1<<" ; err2 = "<<err2;

      //exportResults( u1, u2, cpt);
      exportResults( U1, U2, cpt);

      ++cpt;

    }


  //AUTO( SolAnalytique, cos(M_PI*Px()/2.0)*cos(M_PI*Py()/2.0) );

  //  u1= project( Xh1, elements(mesh1),  SolAnalytique );

  //u2= project( Xh2, elements(mesh2),  SolAnalytique );

  //exportResults( U1, U2, 0);

}


//--//////--//////--//////--//////-------------------------------------------------------------------------------------------------//
//-//--//--//--//--//--//--//--//---------------------------------------------------------------------------------------//
////--//////--//////--//////--////------------------------------------------------------------------------------------------------//

template<uint Dim,uint Order,uint OrderGeo>
void
Localisation<Dim,Order,OrderGeo>::problem1(element_type& U1 , element_type& U2 ) {

  using namespace Life::vf;

  //----------------------------------------------------------------------------------------------------------//

  space_ptrtype Xh=U1.functionSpace();
  mesh_ptrtype mesh=Xh->mesh();

  element_type V(Xh,"v");
  element_0_type v = V.template element<0>();
  element_1_type q = V.template element<1>();

  element_0_type u1 = U1.template element<0>();
  element_1_type p1 = U1.template element<1>();

  element_0_type u2 = U2.template element<0>();
  element_1_type p2 = U2.template element<1>();


  sparse_matrix_ptrtype A( M_backend->newMatrix(Xh,Xh) );
  vector_ptrtype B( M_backend->newVector( Xh ) );
  A1->zero();
  B1->zero();

  //----------------------------------------------------------------------------------------------------------//

  value_type pi = M_PI;
  value_type nu= this->vm()["nu"].template as<value_type>();
  value_type alphaX =  this->vm()["alphaX"].template as<value_type>();
  value_type alphaY =  this->vm()["alphaY"].template as<value_type>();

  value_type penaldir = this->vm()["penaldir"].template as<double>();

  AUTO( Alpha, vec(cst(alphaX),cst(alphaY)) );
  AUTO( grad_Uex , vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) ) );
  AUTO( lap_Uex, -(pi*pi/2.)*cos(pi*Px()/2.)*cos(pi*Py()/2.) );

  //AUTO( f, trans(Alpha)*grad_Uex -nu*lap_Uex);
  AUTO( f,  -nu*lap_Uex);

  //----------------------------------------------------------------------------------------------------------//
    
  form2( Xh, Xh, A, _init=true ) =
    integrate( elements(mesh), _Q<2*Order>(),
	       nu*gradt(u1)*trans(grad(v)) );

  /*  form2( Xh, Xh, A) +=
    integrate( elements(mesh), _Q<2*Order>(),
	       (gradt(u1)*Alpha)*id(v) );
  */

  form2( Xh, Xh, A ) +=
    integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
	       id(v)*idt(p1)
	       +id(q)*idt(u1) );

  form2( Xh, Xh, A ) +=
    integrate( boundaryfaces(mesh), _Q<2*Order>(),
	       - nu*(gradt(u1)*N())*id(v));

  form2( Xh, Xh, A ) +=
    integrate( markedfaces(mesh,mesh->markerName("Bord")), _Q<2*Order>(),
	       //- nu*(grad(v)*N())*idt(u1)
	       + penaldir*id(v)*idt(u1)/hFace());

  form2( Xh, Xh, A ) +=
    integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
	       //- nu*(grad(v)*N())*idt(u1)
	       + penaldir*id(q)*idt(p1)/hFace());


  A->close();


  form1( Xh, B, _init=true ) =
    integrate( elements(mesh), _Q<7>(),
	       f*id(v) );

  form1( Xh, B ) +=
    integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
	       + penaldir*id(q)*idv(p2)/hFace()
	       -id(q)*idv(u2) );



#if 0
  form2( Xh, Xh, A ) +=
    //on( boundaryfaces(mesh) , u, B, cst(0.) ) ;
        on( markedfaces(mesh, mesh->markerName("Bord")) , U1, B, cst(0.) );

  //  form2( Xh, Xh, A ) +=
    //on( boundaryfaces(mesh) , u, B, cst(0.) ) ;
    //    on( markedfaces(mesh, mesh->markerName("Bord")) , p1, B, cst(0.) );
#endif

  this->solve(A,U1,B,false);


}

template<uint Dim,uint Order,uint OrderGeo>
void
Localisation<Dim,Order,OrderGeo>::problem2(element_type& U1 , element_type& U2 ) {

  using namespace Life::vf;

  //----------------------------------------------------------------------------------------------------------//

  space_ptrtype Xh=U1.functionSpace();
  mesh_ptrtype mesh=Xh->mesh();

  element_type V(Xh,"v");
  element_0_type v = V.template element<0>();
  element_1_type q = V.template element<1>();

  element_0_type u1 = U1.template element<0>();
  element_1_type p1 = U1.template element<1>();

  element_0_type u2 = U2.template element<0>();
  element_1_type p2 = U2.template element<1>();

  sparse_matrix_ptrtype A( M_backend->newMatrix(Xh,Xh) );
  vector_ptrtype B( M_backend->newVector( Xh ) );
  A1->zero();
  B1->zero();

  //----------------------------------------------------------------------------------------------------------//

  value_type pi = M_PI;
  value_type nu= this->vm()["nu"].template as<value_type>();
  value_type alphaX =  this->vm()["alphaX"].template as<value_type>();
  value_type alphaY =  this->vm()["alphaY"].template as<value_type>();

  value_type penaldir = this->vm()["penaldir"].template as<double>();

  AUTO( Alpha, vec(cst(alphaX),cst(alphaY)) );
  AUTO( grad_Uex , vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) ) );
  AUTO( lap_Uex, -(pi*pi/2.)*cos(pi*Px()/2.)*cos(pi*Py()/2.) );

  //  AUTO( f, trans(Alpha)*grad_Uex -nu*lap_Uex);
  AUTO( f,  -nu*lap_Uex);

  //----------------------------------------------------------------------------------------------------------//
    
  form2( Xh, Xh, A, _init=true ) =
    integrate( elements(mesh), _Q<2*Order>(),
	       nu*gradt(u1)*trans(grad(v)) );

  /*  form2( Xh, Xh, A) +=
    integrate( elements(mesh), _Q<2*Order>(),
	       (gradt(u1)*Alpha)*id(v) );
  */

  form2( Xh, Xh, A ) +=
    integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
	       -id(v)*idt(p1) );

  form2( Xh, Xh, A ) +=
    integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
	       -id(q)*idt(u1) );

  form2( Xh, Xh, A ) +=
    integrate( boundaryfaces(mesh), _Q<2*Order>(),
	       - nu*(gradt(u1)*N())*id(v));

  form2( Xh, Xh, A ) +=
    integrate( markedfaces(mesh,mesh->markerName("Bord")), _Q<2*Order>(),
	       //- nu*(grad(v)*N())*idt(u1)
	       + penaldir*id(v)*idt(u1)/hFace());

  form2( Xh, Xh, A ) +=
    integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
	       + penaldir*id(q)*idt(p1)/hFace());

  //----------------------------------------------------------------------------------------------------------//

  form1( Xh, B, _init=true ) =
    integrate( elements(mesh), _Q<7>(),
	       f*id(v) );

  form1( Xh, B ) +=
    integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
	       + penaldir*id(q)*idv(p2)/hFace()
	       +id(q)*idv(u2) );

  //----------------------------------------------------------------------------------------------------------//



  A->close();

#if 0
  form2( Xh, Xh, A ) +=
    //on( boundaryfaces(mesh) , u, B, cst(0.) ) ;
        on( markedfaces(mesh, mesh->markerName("Bord")) , U1, B, cst(0.) );

  //  form2( Xh, Xh, A ) +=
    //on( boundaryfaces(mesh) , u, B, cst(0.) ) ;
    //    on( markedfaces(mesh, mesh->markerName("Bord")) , p1, B, cst(0.) );
#endif

  this->solve(A,U1,B,false);

}



template<uint Dim,uint N,uint OrderGeo>
void
Localisation<Dim,N,OrderGeo>::solve( sparse_matrix_ptrtype const& D,
				     element_type& u,
				     vector_ptrtype const& F,
				     bool is_sym )
{

  backend_ptrtype M_backend2(backend_type::build( this->vm() ));
  
  vector_ptrtype U( M_backend2->newVector( u.functionSpace() ) );
  M_backend2->solve( D, D, U, F, false );
  u = *U;
} // Test_Schwarz::solve
 







/**
 * main function: entry point of the program
 */
int
main( int argc, char** argv )
{
  //dim, N, OrderGeo
  Localisation<2,2,1> moncode( argc, argv, makeAbout(), makeOptions() );

  //  moncode.run();
  moncode.runMultiDomain();
  //  moncode.runTest();
}



